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We revisit the problem of determining the real-frequency density response in quantum fluids via 
analytical continuation of imaginary-time quantum Monte Carlo data. We demonstrate that the 
average spectrum method (ASM) is capable of revealing resolved modes in the dynamic structure 
factor of both ort/io-deuterium and liquid para-hydrogen, in agreement with experiments and quan- 
tum mode-coupling theories, while the maximum entropy approach yields only a smooth unimodal 
spectrum. Outstanding issues are discussed. Our work provides the first application of the ASM 
method in realistic off-lattice systems. 
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I. INTRODUCTION 

The computation of real-time correlation functions 
in condensed phase quantum systems is a challeng- 
ing problem for which there is no universally appli- 
cable approach. This fact has necessitated the de- 
velopment of a variety of powerful methods rang- 
ing from exact numerical approaches to approximate 
theories^M&^.a.Vo^biVs^^ie,!? The choice of 

method will in general depend on the type, dimension- 
ality, and size of the system of interest as well as the 
external thermodynamic conditions under consideration. 

The problems that plague exact real-time approaches 
are generally absent in imaginary time.— The dynamical 
sign problem, which renders most direct real time Monte 
Carlo simulations infeasible, is avoided in imaginary time. 
In principle, this means that a well defined Wick rota- 
tion may be used to recover real time information and 
excitation spectral If the precise analytical form of the 
correlation function is known in imaginary time, rotation 
in the complex time plane is straight forward. However 
in most cases of interest the imaginary time correlation 
function must be determined numerically via quantum 
Monte Carlo. In practice, analytic continuation in this 
case is numerically unstable and highly sensitive to sta- 
tistical error.— 

The most commonly used technique for the numeri- 
cal analytical continuation of imaginary time quantum 
Monte Carlo data is the maximum entropy (MaxEnt) 

mc th oc [ 20.21.22.23.24.25.26.27.28.29.30.31.32.33.34.35 J n Ty[ax- 

Ent, the optimal fit of the data is defined in a Bayesian 
manner as the most likely fit that emerges from the com- 
petition between the \ 2 goodness-of-fit and an entropic 
prior S. MaxEnt generally requires an additional means 
to determine the prefactor a ("temperature") of the en- 
tropy term. Once determined, MaxEnt provides a statis- 
tically rigorous and unique fit for the spectral function. 
The strength of the MaxEnt method is that it generally 
provides good estimates of spectral area in the correct 
region of frequency space. However numerous studies 
have shown that, in general, MaxEnt solutions tend to 



be too broad and smooth and often lack clearly defined 
peaks. 3 5JI Other approaches have been devised to over- 
come this "smoothness" problem but these approaches 
have not met with much general success. 

While MaxEnt provides a Bayesian estimate of the 
most likely fit, a different approach may be devised based 
on the notion of averaging over a sequence of possible 
solutions. Such approaches are called stochastic ana- 
lytic continuation methods or average spectrum methods 
(ASM) i 88 ' 39 ' One justification for such a procedure is 
the fact that MaxEnt solutions for different values of a 
may have spurious features that vanish or are diminished 
when averaged. This type of approach was independently 
suggested for the continuation of imaginary time quan- 
tum Monte Carlo data by White and Sandvik ) 38 i 41 while 
similar ideas have been used in other fields . 42 ' 43 While 
these methods have been put forward on rather heuristic 
grounds, the results generated in several studies give sup- 
port to the notion that features not detected by MaxEnt 
may be faithfully revealed. Beach has provided a more 
secure theoretical foundation for such an approach by 
constructing a version of an ASM for which MaxEnt is 
the mean field limit.— Thus, one may view the ASM as 
including fluctuations not incorporated in the MaxEnt 
method. 

In this work, we will use a variant of the ASM ap- 
proach outlined by Slyjuasen 4 ^ to calculate the dynamic 
structure factor of both liquid orf/io-deuterium and liquid 
para-hydrogen. Such cases are useful because extensive 
experimental results exist in these systems that show fea- 
tures that are not captured by MaxEnt. While the ASM 
approach used here is more computationally demanding 
than MaxEnt, we find that this approach is capable of 
revealing such subtle features. Not surprisingly, the ap- 
proach fares better for higher quality of the Monte Carlo 
data. Our work provides the first application of the ASM 
for analytical continuation of quantum Monte Carlo data 
for off-lattice systems with continuous potentials. 

The outline of the paper is as follows: In Sec. IH we 
outline the average spectrum method. Sec. IIIII describes 
the model of liquid ort/io-deuterium and liquid para- 
hydrogen, provides the technical Monte Carlo and in- 
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version details, and presents the results for the dynamics 
structure factor of both liquids. In Sec.|lV]we conclude. 



II. ANALYTIC CONTINUATION AVERAGE 
SPECTRUM METHOD 

The analytic continuation of the intermediate scatter- 
ing functio n 24 ' 33 ' 34 ' 37 ' 45 is based on the Fourier relation 
between the dynamic structure factor S(q,oj) and the in- 
termediate scattering function F(q,t): 
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By performing the replacement t 
detailed balance relation S(q, ~uj) 
tain 



— it, and using the 
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where t, r > 0, and 



F(q, t) = —— Tr 
w ' > ZN 



S(q,w), (2) 
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Here, Z is the partition function, the quantum collective 
density operator is given by p q = X^Li e lq r ° , N is the 
number of particles, q the wave number, and f a is the 
position operator of particle a. 

The reason for introducing the imaginary time inter- 
mediate scattering function, F(q,r), is that, unlike its 
real time counterpart, it is straightforward to obtain it 
using an appropriate path-integral Monte Carlo (PIMC) 
simulation technique ! 46 ' 47 However, in order to obtain 
the dynamic structure factor and the real time interme- 
diate scattering function one has to invert the integral 
in Eq. J2]). Due to the singular nature of the integration 
kernel the inversion of Eq. ([2]) is an ill-posed problem. As 
a consequence, a direct approach to the inversion would 
lead to an uncontrollable amplification of the statistical 
noise in the data for F(q, r), resulting in an infinite num- 
ber of solutions that satisfy Eq. 

In the present study we apply the recently developed 
averaged spectrum method,— related to earlier stochas- 
tic continuation approaches , 38 i 39 i 44 to invert the integral 
in Eq. ([2j) to get S(q, u>). Following the notation of other 
analytic continuation methods, we will refer to F(q,r) 



as the data (input data generated from the PIMC ap- 



proach), K (t, lu) 



as the singular kernel, 



and S(q,cj) as the solution. Furthermore, in what fol- 
lows, we will assume that both the imaginary time axis 
and the frequency axis are discretized such that r = j T Sr 
(j T = 1, • • • N T ) and u = jujSuj (j w = 1, • • ■ N u ). Hence, 
Eq. |(2|) in its discrete form is given by: 



(4) 



where for future reference, the vectors F(q) and S(q) de- 
scribe the data and the solution in discrete space, with 
values given by Fj T (q) and Sj^ (q) , respectively. 

The basic idea behind the ASM is to pick the final so- 
lution for S(q) as the average spectral function obtained 
by averaging over a posterior distribution: 



S(q) = 



fd\S( q )\S( q )P(S(q)\F(q)) 
Jd\S(q)\P(S(q)\F(q)) ' 



(5) 



where S(q) is a solution to Eq. (@]) for a given input 
F(q), and the posterior distribution can be expressed us- 
ing Bayes theorem as i 2Q i 48 



P(S( q )\t(q)) = P(t(q)\S(q))P(S(q)). 



(6) 



In the above, P(S(q)) is the prior probability distribu- 
tion and P(F(q)\S(q)) is the likelihood function. In the 
present study we assume a uniform prior for all positive 
spectral functions with a zero moment that obeys the 



sum rule Ko j 
P(S(q)) as42 



Sj a (q)5u! = Fo(q). Thus, we can write 



P(S(q)) cx 5 \J2K 0Jal S ja ,(q)Sw - F (q) ] J\Q(S ju {q)), 

(7) 

where S(x) is Dirac's delta function and Q(x) is the Heav- 
iside step function. The choice of the prior distributions 
can be tricky. In order to avoid any undesired bias in 
the inversion procedure, we limit the present study to 
the case of a uniform prior that meets the necessary sum 
rule described above. 

The likelihood function appearing in Eq. © describing 
the fluctuations of the imaginary time data is assumed 
to be of a Gaussian form: 



(- n 

P(F(q)\S(q)) cx exp --Tr^F^) - F(g) s W) T E(g)- 1 (F i (g) - F(qf^ 



(8) 



where we have partitioned the PIMC data into n bins each contains m measurements, the vector F l (q) is the 
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Figure 1: Plots of the dynamics structure factor for liquid ort/to-deuterium (left panels) and liquid para- hydrogen (right panels) 
for different values of q. Experimental results are taken from Ref. for liquid ori/io-deuterium and from Ref. [50 for para- 
hydrogen. Values of q are in inverse . 



average result for bin i, and F(q) s ( qS> is the data vector 
corresponding to a specific solution S(g). The width of 
the Gaussian distribution is taken from the covariance 
matrix with elements of given by: 

*w<z) = ^EK(«)-^(«)nw-4(?)). 

i 

(9) 

The above expression for the likelihood function can fur- 
ther be simplified^ to derive a working expression give 
by: 

P(F(q)\S(q)) <xexp{~E(S(q)))} , (10) 
with the "energy" function E(S(q)) given by: 

E(S(q)) = Tr(F( q )~F(q) s ^fi:( q )-\F( q )-F( q ) s M). 

(11) 

For readers who are interested in a more comprehensive 
discussion of the ASM we refer them to Refs. andl-iol. 
Here, we have provided a short overview for completeness 
and note that the averaging procedure described above 
ensures that solutions (S(g)) which over-fit the noise in 
the data (F(q)) are averaged out, and thus the outcome 
is likely to be a smooth spectrum, with realistic features. 

III. MODEL AND RESULTS 

The major goal of the present study is to con- 
duct a direct comparison between the ASM, the 



MaxEnt method, and the QMCT for collective 
density fluctuations in liquid ori/io-deuterium and 
liquid para- hydrogen. Both liquids have been stud- 
ied extensively over the past decade by many different 
techniques H- 30 - 33 - 34 - 37 - 45 - 51 - 52 - 53 - 54 - 55 - 56 - 57 - 58 - 59 - 60 - 61 - 62 - 63 - 64 

These and other studies show that they are characterized 
by quantum dynamical susceptibilities which are not 
reproducible using classical theories. Thus, these liquids 
are ideal to assess the accuracy of methods developed 
for quantum liquids and serve as canonical models in 
the field. 

The static input required by the QMCT and the imag- 
inary time intermediate scattering function required for 
the analytic continuation approaches were generated by 
PIMC simulations at T = 20.7K and p = 0.0254A" 3 
for liquid ori/io-deuterium^ and T = 1AK and p — 
0.0235A" 3 for liquid para-hydrogen.^ 6 - The PIMC simu- 
lations were done using the NVT ensemble with 256 par- 
ticles interacting via the Silvera-Goldman potential ) 67 ' 68 
where the entire molecule is described as a spherical par- 
ticle, so the potential depends only on the radial dis- 
tance between particles. The staging algorithm 6 -^ for 
Monte Carlo chain moves was employed to compute the 
numerically exact input. The imaginary time interval 
was discretized into N T Trotter slices of size e = (3/N T 
with N T = 20 and N T — 50 for liquid orifto-deuterium 
and liquid para- hydrogen, respectively. Approximately 
1 — 2 x 10 6 Monte Carlo passes were made, each pass con- 
sisted of attempting moves in all atoms and all the beads 
that were staged. Data collection was taken every 10 
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Monte Carlo steps, block averaging was done for to = 10 
and the total number of bins was set to n — 1 — 2 x 10 4 . 
The acceptance ratio was approximately 0.25 — 0.3 for 
both liquids. 

In Fig. Q] we show the results for the dynamic struc- 
ture factor S(q,uj) for both liquids and for different val- 
ues of q. The results are compared to the experiments of 
Mukherjee et a/ i 49 ' 70 for ort/io-deuterium and to the ex- 
periments of Bermejo et a/ i 50 ' 56 for liquid para-hydrogen 
for selected values of q (the theoretical values of q are 
slightly different from the experiments due to the limi- 
tations associated with the constant volume simulation 
and finite size effects). Several important features are 
observed experimentally: First, the present of a finite 
frequency peak that disappears around q = 1.4 -1 signi- 
fying the presence of collective density fluctuations (this 
feature is not observed classically for these systems) . Sec- 
ond, the presence of a low frequency peak associated with 
the long-time relaxation of density fluctuations. Finally, 
as q approaches q m ax (qmax ~ 2 _1 is the value of q where 
the static structure factor reaches its first maximum) one 
observes a quantum mechanical de Gennes narrowing of 
the dynamic structure factor. 

These features are semi-quantitatively reproduced by 
the QMCT as shown in Fig. [TJ (a full description of the 
method can be found in our previous publication^ 7 -). In 
particular, the theory captures the position of both the 
low and high intensity peaks for liquid orf/io-deuterium 
and slightly overestimates the position of the high fre- 
quency peak for liquid para- hydrogen. The width of the 
peaks, which are associated with the lifetime of density 
relaxation, are also captured by QMCT as is the quantum 
mechanical de Gennes narrowing of the dynamic struc- 
ture factor. The overall good agreement between the 
QMCT and the experiments is remarkable and signifies 
the importance of the quantum mode-coupling portion 
to the memory kernel at long times . 37 ' 51 ' 53 ' 55 

Turning to the MaxEnt results for S(q,u>) (a full de- 
scription of the calculation can be found elsewhere l37l ) , 
it becomes obvious that MaxEnt fails to provide a quan- 



titative description of the density fluctuations in these 
liquids. It is well known that the MaxEnt approach of- 
ten fails when several distinct spectral features appear 
that are closely spaced in frequency. This is clearly the 
case here where the MaxEnt approach predicts a single 
frequency peak instead of two, at a position that is ap- 
proximately the averaged position of the two experimen- 
tal peaks. Only when the dynamics are characterized by 
a single relaxation time, like the case at g max , does the 
MaxEnt approach provides quantitative results. 

The failure of the MaxEnt result poses a challenge for 
analytic continuation based methods. Can the inversion 
of Eq. |(2]) by other analytic continuation methods pro- 
duce features that are not reproducible by MaxEnt? In 
order to address this problem, we need to average the 
spectrum over the posterior distribution, e.g, solve for 
the average spectrum as given by Eq. ([5]). The averag- 
ing procedure is best done by a Metropolis Monte Carlo 
approach described in Ref. [4y. The basic idea behind 
the Monte Carlo procedure is to start with any guess 
for the spectrum S(q) and then preform a Metropolis 
walk in the solution space to obtain different spectra S(q) 
that are distributed according to posterior distribution 
P(S(q)\F(q)). 

In practice, we follow the recipe of Ref. S3. First, we 
randomly choose a pair of neighboring frequencies, j u 
and j u + 1. We then select a random number ( = [—a, a] 
where a = Sj ul (q)Koj kl + Sj ul ^-i(q)Koj kl ^-i, and make a 
trial move that preserve the sum rule Ko,j u Sj^ {q)Sui — 



S% w (q)=S ju ,(q)+t 
S^q) = S^ +1 (q) - C 



(12) 



The trial move is accepted only if S"™ 1 (q) > with the 
Metropolis probability given by: 



P(S(q) - S n ™(q)) = P(S n ™(q))mm {l,exp [-£ (E(S new (q)) - E(S(q)))] } 



(13) 



where E(S(q)) is given by Eq. (fTTjl . In the present 
study we employ a simple Metropolis algorithm to sam- 
ple the values of S(q) . We did not encounter difficul- 
ties with sampling that require more sophisticated Monte 
Carlo procedures like parallel-tempering^ As a check, 
we started the Monte Carlo procedure from 10 different 
initial conditions. All the different runs converged to the 
same solution within the statistical noise of the Monte 
Carlo sampling. The ASM results shown in Fig. Q] where 
averaged over « 10 s Monte Carlo sweeps, each consist of 



an attempt to change — 512 solution points. 

The main and perhaps, unexpected, result is that the 
ASM captures qualitatively all feature observed exper- 
imentally for the dynamics structure factor, as shown 
in Fig. [TJ Contrary to MaxEnt, the ASM reproduces 
the low and finite frequency peaks, albeit the fact that 
it somewhat underestimates the width of the finite fre- 
quency peak and also slightly overestimates the position 
of the low frequency peak. Overall, the ASM agrees bet- 
ter with the QMCT and with the experimental results 
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Figure 2: Plots of the imaginary time intermediate scattering function for liquid ort/io-deuterium (left panels) and liquid para- 
hydrogen (right panels) for different values of q. Open circles represent the results obtained from the PIMC simulations and 
the solids lines are fits of the ASM to the PIMC results. Values of q are in inverse . 



as compared to MaxEnt. The fact that there are fea- 
tures that are observed experimentally and captured by 
QMCT and ASM and not by MaxEnt method, indicate 
the superiority of the former approaches over MaxEnt. 

Despite the semi-quantitative agreement of the QMCT 
and ASM approaches with experiments, there are still 
significant differences, as clearly is the case for liquid 
para-hydrogen. It still remains a challenge to determine 
the accuracy of both approaches since an exact solution 
for this many-body quantum mechanical problem is still 
(and may always be) out of reach. A direct compar- 
ison with experiments can be misleading, since several 
approximations are made for both QMCT and ASM. 
First, both are based on the assumption that the dy- 
namics evolve over the potential energy taken from the 
work of Silvera and Goldman . 67 ' 68 Furthermore, for sim- 
plicity, para-hydrogen and orifto-deuterium were treated 
as spherical particles, while experimentally this is not 
the case. Finally, for convenience, the simulations were 
done for the NVT ensemble while experiments are con- 
ducted at constant pressure (NPT ensemble). These ap- 
proximation may lead to the difference observed between 
the theory and experiments, but the differences may also 
be attributed to the approximations made by QMCT 
and ASM for the dynamics itself. Nonetheless, the en- 
couraging agreement between QMCT and ASM provides 
stronger support for their accuracy. 

Before we conclude, in Fig. [2] we show the excellent 
agreement between the imaginary data generated by the 
PIMC simulations and imaginary data obtained by the 



inversion of Eq. ([2]). As is well understood, only disagree- 
ment at this level can be used to draw decisive conclusion 
about the inversion process. However, it is quite satisfac- 
tory that the Monte Carlo procedure in the ASM leads 
to an excellent agreement for the imaginary time data. 



IV. CONCLUSIONS 

A comparison between the QMCT, MaxEnt method, 
and ASM has been made for the case of collective density 
fluctuations in liquid ori/io-deuterium and liquid para- 
hydrogen. As far as we know, the present study is the 
first attempt to apply the ASM to dynamical suscep- 
tibilities for off-lattice models. We find that the main 
features observed experimentally for the cases discussed 
above are captured by QMCT and ASM, but not by Max- 
Ent. The results obtained by the ASM appear to coincide 
with the experimental results more closely in the case of 
ori/io-deuterium. In this case the ASM produces results 
that are narrower than those seen in experiment. How- 
ever, it should be expected that the theoretical spectrum 
is indeed sharper than that obtained in experiment be- 
cause no instrumental broadening function is included in 
the theoretical calculations. The better agreement in the 
case of ori/io-deuterium compared to para-hydrogen for 
the ASM may be a result of the better statistical quality 
of the data used in the orf/io-deuterium case. 

Our ASM results for the non-trivial examples pre- 
sented here give strong motivation to continue the in- 
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vestigation of the ASM approach in other systems. One 
important problem worthy of reinvestigation is study of 
density fluctuations in superfluid helium. Boninsegni and 
Ceperley^i studied this problem and found that the sharp 
features associated with roton peaks could not be repro- 
duced by the MaxEnt approach. It would be interest- 
ing to see what the ASM yields for this challenging and 
important problem. Indeed, since the the ASM results 
presented here and by others provide clear evidence that 
analytic continuation based approaches are not limited to 
the case where dynamical susceptibilities are character- 
ized by a single timescale, and since the ASM is a flexible 
approach suitable to the study of any type of correlation 
function for which imaginary time quantum Monte Carlo 
data can be generated, we expect that this approach will 
be further developed and utilized as a powerful means of 
extracting dynamical information from imaginary time 
simulations. 

We would also like to add a word that at this stage it is 
still not entirely clear what should be trusted a priori in 
the ASM results and if the ASM is universally superior 
to MaxEnt. Unpublished work investigating the behav- 



ior of both MaxEnt and the ASM in lattice models of the 
temperature and doping dependence near the Mott tran- 
sition actually suggest that MaxEnt can in some cases 
provide a superior description of certain spectral features, 
but is inferior in its description of others Q For example 
MaxEnt fails completely in describing the Fermi liquid 
behavior of the self energy in the low frequency range 
while, rather surprisingly, the ASM does. Thus, we may 
conservatively suggest that at the very least, the ASM 
provides a complimentary and useful tool for extracting 
real time features from imaginary time quantum Monte 
Carlo data that appears to be superior to MaxEnt at 
least in its ability to uncover sharp spectral features. As- 
sessing the accuracy of the ASM in a more systematic 
manner is another worthy direction of future research. 
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